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Abstract 

The associated production of charged Higgs bosons and top quarks at hadron colliders is an 
^T) I important discovery channel to establish the existence of a non-minimal Higgs sector. Here, we 

-^ ■ present details of a next-to-leading order (NLO) calculation of this process using the Catani- 

^vq . Seymour dipole formalism and describe its implementation in POWHEG, which allows to match 

NLO calculations to parton showers. Numerical predictions are presented using the PYTHIA 

^ ' parton shower and are compared to those obtained previously at fixed order, to a leading order 

H ■ 

Cd ■ calculation matched to the PYTHIA parton shower, and to a different NLO calculation matched to 

the HERWIG parton shower with MC@NLO. We also present numerical predictions and theoretical 

uncertainties for various Two Higgs Doublet Models at the Tevatron and LHC. 
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I. INTRODUCTION 

One of the most important current goals in high-energy physics is the discovery of the 
mechanism of electroweak symmetry breaking. While this can be achieved, as in the Stan- 
dard Model (SM), with a single Higgs doublet field, giving rise to only one physical neutral 
Higgs boson, more complex Higgs sectors are very well possible and in some scenarios even 
necessary. E.g., in the Minimal Supersymmetric SM, which represents one of the most 
promising theories to explain the large hierarchy between the electroweak and gravitational 
scales, a second complex Higgs doublet is required by supersymmetry with the consequence 
that also charged Higgs bosons should exist. 

At hadron colliders, the production mechanism of a charged Higgs boson depends strongly 
on its mass. If it is sufficiently light, it will be dominantly produced in decays of top quarks, 
which are themselves copiously pair produced via the strong interaction. Experimental 
searches in this channel have been performed at the Tevatron by both the CDF [1] and DO 
[2| collaborations and have led to limits on the top-quark branching fraction and charged 
Higgs-boson mass as a function of tan/3, the ratio of the two Higgs vacuum expectation 
values (VEVs), for various Two Higgs Doublet Models (2HDMs). However, if the charged 
Higgs boson is heavier than the top quark, it is dominantly produced in association with top 
quarks with a semi-weak production cross section. The DO collaboration have searched for 
charged Higgs bosons decaying into top and bottom quarks in the mass range from 180 to 
300 GeV and found no candidates |3|]. At the LHC, the ATLAS (and CMS) collaborations 
have already excluded top-quark branching ratios to charged Higgs bosons with masses of 90 
(80) to 160 (140) GeV and bottom quarks above 0.03-0.10 (0.25-0.28) using 1.03 fb'^ (36 
pb^^) of data taken at y/S = 7 TeV 4,|5|. At 14 TeV and with an integrated luminosity of 30 
fb~\ the discovery reach may be extended to masses of about 600 GeV using also the decay 
into tau leptons and neutrinos J6l, 7l|. It may then also become possible to determine the spin 
and couplings of the charged Higgs boson, thereby identifying the type of the 2HDM realized 
in Nature. Searches for pair-produced charged Higgs bosons decaying into tau leptons and 
neutrinos, second generation quarks and VT-bosons and light neutral Higgs bosons have been 
performed at LEP and have led to mass limits of rriH > 76.7 (78.6) GeV for all values of 
tan/3 in Type-I (8| (Type-II [91]) 2HDMs, where only one (both) Higgs doublet (s) couple to 
the SM fermions. Indirect constraints from flavor- changing neutral currents (FCNCs) such 



as 6 — )■ S7 can be considerably stronger, e.g. tuh > 295 GeV for tan /3 > 2 in the absence of 



other new physics sources 10|. 



In this paper, we concentrate on the associated production of top quarks and charged 
Higgs bosons at hadron colhders, which is of particular phenomenological importance for a 
wide range of masses and models. Conversely, s-channel single, pair, and associated pro- 
duction of charged Higgs bosons with VT-bosons are less favorable in most models. Isolation 
of this signal within large SM backgrounds, e.g. from top-quark pair and W^-boson asso- 
ciated production, and an accurate determination of the model parameters require precise 
predictions that go beyond the next-to-leading order (NLO) accuracy in perturbative QCD 



obtained previously ll|-|l3| . We therefore present here details of our re-calculation of this 



process at NLO using the Catani-Seymour (CS) dipole subtraction formalism |14j. Il5|. The 
virtual loop and unsubtracted real emission corrections are then matched with a parton 
shower (PS) valid to all orders in the soft-coUinear region using the POWHEG method 
|l6l . llTj in the POWHEG BOX framework 18|. A similar calculation has been presented 
using the Frixione-Kunszt-Signer (FKS) subtraction formalism [19] and matching to the 



HERWIG PS with the MC@NLO method 



2o|. 



Other new physics processes recently imple- 



mented in MC@NLO include, e.g., the hadropro duct ion of additional neutral gauge bosons 



2l| . Unlike MC@NLO, POWHEG produces events with positive weight, which is important 
when the experimental analysis is performed via trained multivariate techniques. POWHEG 
can be easily interfaced to both HERWIG [22] and PYTHIA [23| and thus does not depend 
on the Monte Carlo (MC) program used for subsequent showering. 

The remainder of this paper is organized as follows: In Sec. HIl we present details of 
our NLO calculation of top-quark and charged Higgs-boson production. We emphasize 
the renormalization of wave functions, masses, and couplings, in particular the one of the 
bottom Yukawa coupling, in the virtual loop amplitudes as well as the isolation and can- 
cellation of soft and collinear divergences with the Catani-Seymour dipole formalism in the 
real emission amplitudes. The implementation in POWHEG is described in Sec. Illli As 
the associated production of top quarks with charged Higgs bosons is very similar to the 



one with PF-bosons 2J], we concentrate here on the differences of the two channels. We 
also emphasize the non-trivial separation of the associated production from top-quark pair 
production with subsequent top-quark decay in scenarios, where the charged Higgs boson is 
lighter than the top quark, using three methods: removing completely doubly-resonant dia- 
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FIG. 1: Tree- level diagrams for the associated production of charged Higgs bosons and top quarks 
at hadron colliders in the s-channel S and the t-channel T. 

grains, subtracting them locally in phase space, or including everything, so that top-quark 
pair production with the subsequent decay of an on-shell top quark into a charged Higgs 
boson is effectively included at leading order (LO). In Sec. llVt we present a detailed numer- 
ical comparison of the new POWHEG implementation to the pure NLO calculation without 

n n 

PS [12!], to a tree-level calculation matched to the PYTHIA parton shower [25j; and to the 



20[. We also give numerical predictions 



MC@NLO implementation with the HERWIG PS 

and theoretical uncertainties for various 2HDMs at the Tevatron and LHC Our conclusions 

are summarized in Sec. |Vl 



II. NLO CALCULATION 



A. Organization of the calculation 



At the tree level and in the five-fiavor scheme with active bottom (6) quarks as well as 
gluons {g) in protons and antiprotons, the production of charged Higgs bosons {H~) in 
association with top quarks (t) occurs at hadron colliders via the process h{pi) + g{p2) -^ 
H~{ki) + t{k2) through the s- and t-channel diagrams shown in Fig. [1] The massive top 
quark is represented by a double line, whereas the bottom quark is treated as massless and 
represented by a single line. The Born matrix elements can then be given in terms of the 
Mandelstam variables 



S = {Pl+P2f = {h + k2f, 

t = (p2 - k2y = {ki - pif, and 



u = ml + m^ — s — t. 



(1) 
(2) 
(3) 



A NLO calculation in a four-flavor scheme, where the bottom quark is treated as massive 
and generated by the splitting of an initial gluon, has been presented elsewhere [26|, but 
the effect of the bottom mass through the parton densities was subsequently found to be 
strongly suppressed compared to its impact on the bottom Yukawa coupling [27]. The four- 
momenta of the participating particles have been ordered in accordance with the POWHEG 
scheme, where the initial-state particles with four-momenta pi and p2 are followed by the 
four-momentum ki of the final-state massive colorless particle and then the four-momentum 
^2 of the outgoing massive colored particle. The additional radiation of a massless particle 
(gluon or light quark), that occurs at NLO in real emission diagrams, is assigned the last 
four- momentum k^. 

The hadronic cross section 



(^AB 



y2 / dXa fa/A{Xa, f^p) / dXb fb/B{Xb, fJ-p) (^{Pl, P2] fJ'p) 

„ I, Jo Jo 



(4) 



is obtained as usual as a convolution of the parton density functions (PDFs) fa/A,b/B{xa,b, /^f) 
with the partonic cross section 

^(pi,P2;/u|) = fT^''(pi,P2) +cr'^''''(pi,P2;/i|) (5) 

with partonic center-of-mass energy s = XaX^S, S being the hadronic center-of-mass energy. 



Its LO contribution 



^^^(Pl,P2) 



2s 



d^^^^lM 



Born 



(6) 



is obtained from the spin- and color- averaged squared Born matrix elements |A^bc 
through integration over the two-particle phase space 6?$*^^-' and flux normalization. 



B. Virtual corrections and renormalization scheme 



Like the LO partonic cross section a^^, its NLO correction 
^^"^(pi,P2;/.|) = a^^«^^}(pi,p2) + a^^^^^>(pi,p2) 



+ dx 
Jo 

has a two-body final-state contribution 



a^^^{2>(x; xp^,p,; /i|) + a^^^^^>(x;pi, xp,; /i|) 



a 



NLO{2} 



{Pl,P2) 



[rfa^(pi,P2) +rfa^^(pi,P2) ®I],^0 



rf<l>(2) 



2 Re 



A^i-ioopA^Born + 2{H,t;b,g I 1(e) I H,t;b,g)2 



e=0 



(7) 



(8) 



which consists of the virtual cross section da^ , i.e. the spin- and color-averaged interference 
of the Born diagrams with their one-loop corrections, and the Born cross section da^'^ 
convolved with a subtraction term I, which can be written as 2{H., t;b,g \ 1(e) | H, t; b, (7)2 and 
which removes the infrared singularities present in the virtual corrections. The three-body 
final-state contribution a^^'~'^^^ and the finite remainders a^^'^^^^{x, ...) of the initial-state 
singular terms will be described in the third part of this section. 

The ultraviolet divergencies contained in the virtual cross section da^ have been made 
explicit using dimensional regularization with D = 4 — 2e dimensions and are canceled 
against counterterms originating from multiplicative renormalization of the parameters in 
the Lagrangian. In particular, the wave functions for the external gluons, bottom and top 
quarks are renormalized in the MS scheme with 



SZ, 



2N, 



c 



SZ, 



b,t 



as 



-—Cp/^uv, 
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-N, 



c 



-N, 



A 



uv 



and 



(9) 
(10) 



where A[/y = 1/e — 7^; + ln47r, •je is the Euler constant, Nq = 3 and A''^;' = 6 are the 
total numbers of colors and quark flavors, respectively, and Cp = {N^, — l)/{2Nc). The 
counterterm for the strong coupling constant as = g\/{A'K) 



9s 



iTC 



Auvi^^c 



-Nf 
3 



?ln4 

3 ml 



(11) 



is computed in the MS scheme using massless quarks, but decoupling explicitly the heavy 
top quark with mass rrit from the running of as [28]. The top-quark mass entering in the 
kinematics and propagators is renormalized in the on-shell scheme, 

asifi%).^ f . 4 "2 



rat 47r V 3 m^/ 



(12) 



On the other hand, we perform the renormalization of both the bottom and top Yukawa 
couplings in the MS scheme. 



Jb,t 



''R) 



3CfA 



uv- 



(13) 



This enables us to factorize the charged Higgs-boson coupling at LO and NLO, making the 
QCD correction (K) factors independent of the 2HDM and value of tan/3 under study. In 
particular, in Eq. ( fT3l) we do not subtract the mass logarithm, but rather resum it using the 



running quark masses 



mQifiR) = mQ{MQ)- 



c[as(fiR)/n 



c[asiMQ)/TT 



(14) 



in the Yukawa couplings, where 



/23 \ 12/23 

c{x) = (—x) (1 + 1.175X + l.SOlx^ 



(15) 



for rub < fJ.R < rrit and 



'7 \4/7 



c{x) = i-xj (1 + 1.398X + 1.793x^) 



(16) 



for fiji > rub^t 29|. The starting values of the MS masses are obtained from the on-shell 
masses Mq through the relation 

Mo 



rfiQiMQ) 



1 + 



4"s(A£q) 



■K V "^ / 



(17) 



with Kb ^ 12.4 and K, ^ 10.9 



30,131 



After the renormalization of the ultraviolet singularities has been performed as described 
above, the virtual cross section contains only infrared poles. These are removed with the 
second term in Eq. ([8]), i.e. by convolving the Born cross section with the subtraction term 



Q^ 



1(e) = I2{e,^^]{k2,mt]) + lb{e,^i^]{k2,mt],pi) + lg{e, ^^]{k2,mt],p2) + Ibg{e, fi^;pi,p2), 

(18) 
where in our case l2(e,/i^; {k2,mt}) = 0, since there are no QCD dipoles with a final state 
emitter and a final state spectator. The dipoles depending on one initial-state parton (a = 
b,g) with four- momentum pi [i = 1,2) are 

as {AttY f 1_ 

-12 
■t 



I„(e,/i ■,{k2,mt},Pi) 



27rrri-e) 1 T?'^*''^" 



T?f£) (Vt{s,a,m„0;e)-^ 



+ rT-i2 -'-a ■ -'-* 



Li 

+ rt(/i, mt; e) + 7i In h 7t + Kt 

Sta 



Sat 



(19) 



where Ta,t denotes the color matrix associated to the emission of a gluon from the parton a or 
the top quark t, the dimensional regularization scale /i is identified with the renormalization 
scale fj,R, and Sta = Sat = '^Pik2- The kernels 

Vg{sgt,0,mt;e,2/3) = V^''\sgt,0,mt;e) + Vf''\sgt,0,mt;2/3) 
consist of the singular terms 



(20) 

(21) 
(22) 



m} 1 , n 771? TT 

^ - - In'' — 

Sta 4 Sta 

With Ql = Q 



A + -ln^--ln^^----ln^ln|^-iln^ln|^ (23) 
2e2 2e s,„ 4 s,„ 12 2 s,, Qi 2 gf„ Qf, ^ ^ 



-ai 



,{Af5) 



Sta + f^t + '^a ^^d the non-singular terms 

2 



Vr'^(5ta,mt,0) 



-l^T^ + y 



T,^ Qta 



Lig 



Sta 
Qla 



-2 In 



Sfa 



Q^ 



2 2 

— ^In 



ta 



V, 



^^^^(Sbt,0,mt) = ;j^ 



In 



Qlt 



21n^^~"^*-2 



"^t 



Q 



fet 



Vf^)(V,0,mt;2/3) = ^ 

9 



T'2 
S 

SAT, 



lni^-21n^^i^^-2 



Q^ 



Qgt -mt nit 



Qbt + mt_ 
mt 

Qgt+ITlt 



Sta Q 

^2 



2 ' 
ta 



"6 



Sbt 



(24) 



+ --Li,(^),(25) 



"" T i ( '^* 



+ -T^ In -^'^^ '^ + 7^-^ o- (26) 

Qgt Qgt + mt 3 J 

The constant k in Eq. f lT9|) is a free parameter, which distributes non-singular contributions 
between the different terms in Eq. ([7]). The choice n = 2/3 considerably simplifies the gluon 
kernel. For massive quarks, one has in addition 



Tt{fi,mt;e) = Cp 



1 , m? 

- n ^ - 2 

2 /x2 



(27) 



while 



7g 



2^^ 



73 



11 2 

— iVc - -TnNf 
6 3 



(28) 



and 



K„ 



~6 



Cf 



K„ 



67 
18 






A^r. - —Tj,N 



^c 



9 



iJ^v/ 



(29) 



with Tr = 1/2 and Nj = 5 the number of light quark flavors. The last term in Eq. ( IT8 



Ifeg(e,/i^;pi,P2) 



a.o 



(4vr) 



27rr 1-e 



rp2 ■'■9'T'' 
9 



+ {g^h) 



f^^ 



Sbg 



Ig 



IT 



7- T# + ^ -T^ + 7. + ^. 



(30) 



depends on both initial-state partons. Since the process we are interested in involves only 
three colored particles at the Born level, the color algebra can be performed in closed form. 
To be concrete, we have 

T,-T,\H,t-b,gh = -(cF-^yH,t-b,gh = ^\H,t-b,gh, (31) 

Nn 

T,,fTg\H,t;b,gh = -^\H,t;b,g)2, (32) 

Tlt\H,t;b,g)2 = Cp \H,t; b, g)2, e^nd (33) 

Tl\H,t-b,g)2 = Nc\H,t;b,g)2. (34) 



C. Real corrections 

The second term in Eq. (^^ 
a''''''^''\puP2) = J d<!>^'^[\M3,i,{h,k2,ks;puP2W- Yl 1^ih,k2,h;puP2)} (35) 

dipoles 

includes the spin- and color-averaged squared real emission matrix elements 



|-^3,ij(^i; ^2, ^3;Pi,P2)p with thrce-particle final states and the corresponding uninte- 
grated QCD dipoles V, which compensate the integrated dipoles I in the previous section. 
Both terms are integrated numerically over the three-particle differential phase space d^^^\ 
The real emission processes can be grouped into the four classes 

(a) b{pi) + g{p2) -^ H-{h) + t{k2) + g{k3), 

(b) g{pi) + g{p2) ^ H-{k^) + t{k2) + 6(^3), 

(c) b{pi) + q/q{p2) -^ H-{ki) + t{k2) + q/q{k^), and 

(d) g(pi) + q{p2) ^ H-{k,) + t{k2) + 6(A;3), 

where the second process (b) can be obtained from the first one (a) by crossing the four- 
momenta ^3 and —pi and multiplying the squared matrix element by a factor of (—1) to take 
into account the crossing of a fermion line. The processes in the two other classes (c) and 
(d) can interfere when q = b, but these contributions are numerically negligible due to the 
comparatively small bottom quark parton distribution function. Process (d) is furthermore 
convergent for q = u,d,s and c. 



The sum over the dipoles in Eq. (135!) includes initial-state emitters ab with both initial- 
and final-state spectators c (D"*'^ and I^"*) and the final-state emitter ab with initial-state 
spectators c (^Jjf,). For the three divergent processes, we have 

(a) : J2 ='D^''' + ^3''^ + ^1' + Vf + V\g + V%, (36) 

dipoles 

(b) : Y^ = p9iM2 + p<?2Mi + pf ^ + Pf ^ and (37) 

dipoles 

(c) : ^ =2)9«.fc + pf. (38) 

dipoles 

Denoting by a the original parton before emission, b the spectator, and i the emitted particle, 
the dipole for initial-state emitters and initial-state spectators is given by 

^P"*'" = -1^^ 2 AH. i- ai, b I ZA-^V'^*'" I H, i; at, b)2,ab, (39) 

where the momentum of the intermediate initial-state parton ai is p^^ = Xi,ab Pa with Xi^ab = 
{PaPb — kiPa — kiPb) / (PaPb) , the momcutum pb is unchanged, and the final-state momenta kj 
with j = 1, 2 are shifted to 

with K^ = Pa ^ Pb ~ k'i and K^ = p'^- + p'^. The necessary splitting functions V^*'^ for 
{ai,b} = {qg, g', gg.q] gq, g',qq,q} can be found in Ref. [ij]. The dipole for initial-state 
emitters and a final-state spectator, which is in our case the top quark t, is given by 

^" = -^V- 2,a.{Hrt-ai,b\ ^^4f^Vr I H,i;m,b),^-,, (41) 

where the momentum of the intermediate initial-state parton ai is p^^ = x^^aPa with xa^a = 
[Paki + PaPt — kiPt)/{paki + PaPt), the momcutum pb is unchanged, and the momentum of 
the final-state top quark pt is shifted to p^ = k^ +p^ — {1 — Xit^a)Pa- The necessary splitting 



functions Vf* for {ai,t} = {qg,t; gg,t; gq,t;qq,t} can be found in Ref. [l5[. Finally, the 



dipole for final-state emitter (the top quark t) and initial-state spectator a is given by 

K = -tAt^ ^AH, tt- a, b I ^-^Yl I H, h- a, 6)2,a, (42) 

where the momentum of the initial parton a is shifted to p^ = x^^aPa with x^^a = 
{PaK + PaPt — kiPt)/{paki + PaPt), the momcutum p^ is unchanged, and the momentum 

10 



of the intermediate final-state top quark pt is pf^ = fcf + pf — (1 — Xit^a)Pa- The required 
sphtting function V^^ can again be found in Ref. [l5| . 

The last terms in Eq. (jTj) are finite remainders from the cancellation of the e-poles of the 
initial-state collinear counterterms. Their general expressions read 

/ dxa''''''^^^x;xpi,p2;ixl) = V /" dx / fe,^ (a;pi,p2) ® (K + P)"'"' (x)] (43) 

Jo ^ Jo J2^ ^ ^=0 

= y^ dx d^^^\xpuP2) 2,a'b{ki,k2;xpi,p2\K''^'''{x) + P''''''{x;i/p)\ki,k2;xpuP2)2,a'b 
, Jo J 



and similarly 
given in Ref. 



or (a -H- b) and (pi ^ ^2)- The color-charge operators K and P are explicitly 



15|- 



III. POWHEG IMPLEMENTATION 



The calculation in the previous section has been performed using the Catani-Seymour 
dipole formalism for massive partons |14j.ll5|. For the implementation of our NLO calculation 
in the POWHEG Monte Carlo program, we need to retain only the Born process, the finite 
terms of the virtual contributions, and the real emission parts of our calculation, since all 
necessary soft and collinear counterterms and finite remnants are calculated automatically 
by the POWHEG BOX in the FKS scheme [l9|. Soft and collinear radiation is then added 
to all orders using the Sudakov form factor. In this section, we briefly describe the three 
relevant contributions, following closely the presentation in Ref. [I81], and address the non- 
trivial separation of the associated production of charged Higgs bosons and top quarks from 
top-quark pair production with subsequent top-quark decay in scenarios, where the charged 
Higgs boson is lighter than the top quark. 

A. Born process 

In the POWHEG formalism, a process is defined by its particle content. Each particle 



is encoded via the Particle Data Group numbering scheme [3^ except for gluons, which are 
assigned the value zero. The order of the final state particles has to be respected. Colorless 
particles are listed first, then heavy colored particles, and finally massless colored particles. 
The Born processes (two with respect to the different bg and gb initial states) are defined 
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TABLE I: Integration limits for the hadronic cross section. 



Variable V 


Vmin 


y-max 


T 

y 
t 


S 

ilnr 
\{ti-t2) 


1 

-ilnr 

i(tl+t2) 










ti = rrif + rnjj — s, t2 


= V(^- 


- "it - mj^y - 


- Amfm'j^ 



with f lst_nborn = 2 and are listed as 



{hg^ R-t) = [5,0,-37,6] 



and 



(44) 



(45) 



[gh^ R-t) = [0,5,-37,6] 

in the subroutine init_processes. 

In the subroutine born_phsp, the integration variables xborn(i) for the Born phase space 
are generated between zero and one. The hadronic cross section is then obtained from the 
differential partonic cross section da via the integration (see Eq. (j4j)) 



O-AB = X^ / dXafa/A / dXbfb/B / -Trdt 

dr / dy fa/A fb/B / 



dt 



dt. 



(46) 



where ft/i is the PDF of parton i inside hadron / with momentum fraction Xi and where we 
have performed the change of variables 



In 



Xn 



and r = XaXb- (47) 

y/XaXb 

The integration limits are given in Tab.[Tl The Jacobian for the change of integration variables 
from xborn(i) to (r, t/,t) 



A 



jac 



[T„ 



'Tmin) ^ yVmax ymin) ^ l,^ 



max ^mm 



Zrt 



(48) 



has to be multiplied with 27r for the integration over the azimuthal angle 0, which is randomly 
generated by POWHEG. The different kinematical variables can then be constructed in the 
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FIG. 2: Color flow in the Born contribution [5,0,-37,6] and for switched incoming partons 
[0,5,-37,6]. 

center-of-mass reference frame as well as in the laboratory frame via boosts. The renormal- 
ization scale fiR and factorization scale ^p are set in the subroutine set_f ac_ren_scales 
according to the usual convention 

flR = IJ^F = -^ , (49) 

where k is to be varied around two for uncertainty studies. Both the born_phsp and the 
set_f ac_ren_scales subroutines can be found in the file Born_phsp.f. 

All other routines relevant to the Born process are contained in the file Born.f. The 
subroutine setborn contains the factors for the color-correlated Born amplitudes, which are 
related to the Born process through the color factors quoted in Eqs. (13T]) and (l32l) . The 
subroutine borncolor_lh contains the color flow of the Born term in the large- A'"c limit 
shown in Fig. [2j The routine compborn contains the spin-correlated Born matrix element 



MZrn 



- (s^s'' + 5^r^ + r^5^ + m''\ (so) 

before summing over the initial gluon polarizations as well as 

A^Born = -9f.uMZrn, (51) 

where Qf^i, is the metric tensor. 

B. Virtual loop corrections 

The renormalized virtual cross section is defined in dimensional regularization and in the 
POWHEG convention by 



V 






13 



(52) 



TABLE II: Process numbers of the different real emissions. Here q = d,u, s, c. 



Process number Initial state 


Process number 


Initial state 


1 bg 


16-19 


qb 


2 gb 


20-23 


QQ 


3 gg 


24-27 


QQ 


4-7 bq 


28 


bb 


8-11 qb 


29 


bb 


12-15 bq 


30 


bb 



where |A^BornP is now the squared Born matrix element computed in D = 4 — 2e dimensions 
and where the remaining double and simple infrared poles are proportional to 



C2 



1 3^^ 

-—- JMc and 

2Nc 2 ^ 



5 - 4 In * 



mf 



Nc 
12 



-37 



121n^ + 121n^^^)+iiV^. 
mf mf / 6 



(53) 
(54) 



ANc ^ llOf / ±^ X 110^ IILf 

The POWHEG implementation needs then only the finite coefficient V/j„, which has been 
organized into terms stemming from scalar 2-, 3- and 4-point integral functions Bq, Cq and 
Do plus remaining terms and can be found in the file virtual . f . Non-divergent Co-functions 
and Euler dilogarithms are computed using routines contained in the file loopfun.f. 



C. Real emission corrections 

In the subroutine init_processes, the index of the first colored light parton in the 
final state is defined, which is in our case the additional jet from the real emission 
(f lst_lightpart = 5). All flst_nreal = 30 real emission processes are then assigned 
a number according to the list given in Tab. [TTl The expressions of the squared real emission 
matrix elements are given in the file real_ampsq.f . 

D. Separation of associated production and pair production of top quarks 

If the charged Higgs-boson mass mn is lower than the top-quark mass rrit, the antitop 
propagator of the real emission amplitudes shown in Fig. |3] can go on shell, resulting in a 
drastic increase of the total cross section. In other words, the prevalent production mecha- 
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T3G0CDOT?< 



FIG. 3: Real emission contributions in the gluon-gluon and quark-antiquark channels with an 
antitop propagator that can go on shell. 

nism becomes the on-shell production of a ti pair, followed by the decay of the antitop quark 
into a charged Higgs boson. The corresponding Feynman graphs contribute to top-antitop 
production at LO with the charged Higgs boson being produced in top-quark decays, but 
also to tH~ production at NLO. The relevant NLO processes are free from coUinear and soft 
singularities. 

At this point the problem arises how to separate the two production mechanisms. In 
the literature two methods have been proposed: Diagram Removal (DR) and Diagram Sub- 
traction (DS) 33|. Both remove the top-quark resonance from the cross section, but the 
procedure for combining top pair production with the associated production is not com- 
pletely clear. If we separate the amplitudes of a real emission process with colliding partons 
a and b into contributions J^^b^ which proceed through tt-pro duct ion, and contributions 
A^^f^ , which do not. 



4tH- 



squaring the amplitudes gives rise to three different quantities: 



(55) 



\M 



ab\ 



\MT\' + 2Re{MTM'::) + \M':, 



tt 1 2 



Sab + '^ab + T^. 



■'ab- 



(56) 



The term Vab contains neither collinear nor soft singularities, while the interference term lab 
contains integrable infrared singularities. These terms are therefore sometimes referred to 
as subleading with respect to those in Sab, which contains all infrared singularities and must 
be regularized, e.g., via the subtraction formalism. DR requires removing ti production 
at the amplitude level. The only contributing element is then Sab- Since it contains all 
divergencies, the dipoles used in the rriH > nit case remain valid. In the DS scheme, one 
subtracts from the cross section the quantity 



da' 



sub ^ fBwjmH- 



J[(ti) 



(57) 
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locally in phase space. The momenta are reorganized so as to put the i quark on its mass 
shell. Although gauge invariant, this procedure is still somewhat arbitrary. We therefore 
introduce here a third option, where nothing is removed or subtracted from the associated 
production, but simply the full production cross section is retained. Once a sample of events 
is generated, one can then still decide to remove events near the resonance region and replace 
them with events obtained, for example, with a full NLO implementation of ti production. 

In our POWHEG code, we implemented all three methods described above. DR is the 
simplest case. If the flag DR is set to one in the file powheg. input, the resonant diagrams of 
Fig. [3] are simply not included. For the other two procedures, i.e. DS and keeping the full 
cross section, DR should be set to zero. The s-channel propagators of the i quark in the real 
amplitudes are then replaced by a Breit-Wigner form. Setting the flag DS to one turns on 
diagram subtraction. If neither DS nor DR are set to one, the full cross section is computed. 
In this case it is, however, hard to probe the i pole with sufficient accuracy in the Monte 
Carlo integration. An additional flag sepresonant is therefore introduced that, when set 
to one, causes POWHEG to treat the resonant contributions as a regular remnant. This is 
possible since they do not require subtractions. A specific routine for the generation of the 
phase space of the regular remnant ensures that appropriate importance sampling is used 
in the i resonant region. 

While with the DR or the full scheme the fraction of negative weights is very small, this 
is not the case in the DS scheme. Here the real cross section can become negative in certain 
kinematic regions, so that POWHEG must then be run with the fiag withnegweights set 
to one. Negatively weighted events are then kept, but are hard to interpret, since they 
correspond to the subtraction of an ad hoc quantity from the cross section. 

Removing diagrams at the amplitude level causes the loss of gauge invariance. A con- 



siderable part of Ref. [33| has been dedicated to the analysis of the corresponding impact 
on Wt production. There, different gauges were considered for the gluon propagator, and 
differences at the per-mille level were found. Note, however, that gauge invariance is not 
only spoiled through the gluon propagator, but also when the polarization sum 

P^''{k)= J2^'ik,^y{k,X) (58) 

A=l,2 

of external gluons is replaced by 

P^"^(/t) = -g^'•' (59) 
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for simplicity. Here, k^^ is the four-niomentum, A is the polarization, and e^{k,X) is the 
polarization vector of the external gluon. Eq. fl59|) includes not only physical transverse, but 
also non-physical gluon polarizations that must be canceled by ghost contributions. Remov- 
ing individual diagrams then causes the loss of gauge invariance. We therefore abandon the 
use of the simple polarization sum, Eq. f l59p . and sum instead only over physical states with 

P^-'^k) = -g^" ^ [r]^k^'k'' -k-T] {k^'f]'' + r]^'k'')] , (60) 

where ■i]'^ is an arbitrary four- vector transverse to the polarization vector e^. When calcu- 
lating a gauge invariant quantity, the ?7-dependence would drop out, but this will not be 
the case in DR as argued above. For the channels with two external gluons and incoming 
four- momenta pi and p2, we choose for the polarization vectors 

Vi = P2 and 1J2 = pi- (61) 

IV. NUMERICAL RESULTS 
A. QCD input 

ladrons, we use the set CTIO 



For the parton density functions (PDFs) in the external 
obtained in the latest global fit by the CTEQ collaboration 3J]. It has been performed at 
NLO in the MS factorization scheme with nf = 5 active flavors as required by our calculation. 
The employed value of as{Mz) = 0.118, close to the world average, is equivalent to setting 
the QCD scale parameter A^^ to 226.2 MeV as in the previous flts. We also adopt their 
value for the bottom quark mass of mf, = 4.75 GeV and for the top-quark mass of rrit = 172 
GeV and not the newest average value oirrit = 173.2 GeV obtained in direct top observation 
at the Tevatron 35|], as the former value corresponds nicely to the one in the MC@NLO 
publication 20| that we will compare with later in this section. For the sake of easier 
comparisons we also adopt the default scale choice fip = A^i? = {f^H + ^7ii)/2 as in the 
MC@NLO study. We use the versions HERWIG 6.5.10 and PYTHIA 6.4.21 with stable 
top quarks and Higgs bosons and no kinematic cuts to simplify the analysis. Multiparticle 
interactions were neglected. For a discussion of the numerical impact of the bottom mass 



in the PDFs we refer the reader to Ref. 27 1 
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B. Two Higgs Doublet Models 

New particles with masses in the TeV range, that couple to quarks at the tree level, can 
strongly modify the predictions for Flavor Changing Neutral Current (FCNC) processes, 
since these are absent at tree level in the SM. Thus all extensions of the SM, including 
the 2HDMs, must avoid conflicts with the strict limits on FCNCs, such as the electroweak 
precision observable Rb = T{Z — t- hb)/T{Z — > hadrons) or the branching ratio BR(i? — >■ 
Xs'-)). In 2HDMs, tree- level FCNCs are traditionally avoided by imposing the hypothesis 
of Natural Flavor Conservation (NFC), which allows only one Higgs field to couple to a 
given quark species due to the presence of a flavor-blind Peccei-Quinn f/(l) symmetry or its 



discrete subgroup Z2 36|. Alternatively, all flavor- violating couplings can be linked to the 
known structure of Yukawa couplings and thus the Cabibbo-Kobayashi-Maskawa (CKM) 
matrix under the hypothesis of Minimal Flavor Violation (MFV) [37]. Both hypotheses 
have recently been compared with the result that the latter appears to be more stable under 



quantum corrections, but that the two hypotheses are largely equivalent at tree level [38 1. 
In a general 2HDM, one introduces two complex SU(2)-doublet scalar fields 

( i>t \ 

$i = with 2 = 1,2, (62) 

VK+0r+^#)/V2; 

where the Vacuum Expectation Values (VEVs) f 1^2 of the two doublets are constrained by 
the ly-boson mass through t>^ = t>^ + t>| = Am^/g'^ = (246 GeV)^ [39]. The physical charged 
Higgs bosons are superpositions of the charged degrees of freedom of the two doublets, 

H^ = -sin/3 0^ + cos/3 0^, (63) 

and the tangent of the mixing angle tan/3 = ^2/^1, determined by the ratio of the two 
VEVs, is a free parameter of the model, along with the mass of the charged Higgs bosons 
rriH- The allowed range of tan/3 can be constrained by the perturbativity of the bottom- 
and top-quark Yukawa couplings {yt^ < 1) to 1 < tan/3 < 41. Note that in the Minimal 
Supersymmetric SM (MSSM) rriH > myy at tree level. The possible assignments of the 
Higgs doublet couplings to charged leptons, up- and down-type quarks satisfying NFC are 
summarized in Tab. Illli 

In the Type-I 2HDM, only $2 couples to the fermions in exactly the same way as in the 
minimal Higgs model, while $1 couples to the weak gauge bosons |40|. The Feynman rules 
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TABLE III: Couplings of the two Higgs doublets $1^2 to up-type quarks (n), down-type quarks 



(d), and charged leptons (/) in 2HDMs satisfying Natural Flavor Conservation 



4ii. 



Model 



$1 
$2 



Type-I Lepton-specific Type-II Flipped 



n,(i, / 



I 



d,l 



u 



d 
uA 



for the charged Higgs-boson couphngs to quarks in this model, with all particles incoming, 



are 



H'^Uidj : 



W 



-Vij{coi /3 rriu.PL - cot /3 nxd^PR) 



(64) 



where Vij is the CKM matrix and Pl,r = ( 1 =F 75 )/\/2 project out left- and right handed quark 
eigenstates. As can be seen from Tab. IIIIl these couplings are the same in the lepton-specific 
2HDM. 

In the Type-II 2HDM, $2 couples to up-type quarks and $1 to down-type quarks and 
charged leptons. The Feynman rules for charged Higgs-boson couplings to quarks in this 
model, with all particles incoming, are 



W 



H^Uidj : ^ 



- Vij (cot /3 rriu, Pl + tan /3 m^^, Pr) . 



(65) 



As can again be seen from Tab. IIIIl they are identical to those in the flipped 2HDM 41 1. 

Since the NFC and MFV hypotheses allow for the possibility that the two Higgs doublets 
couple to quarks with arbitrary coefficients A\ ^, there exists also the possibility of a more 
general 2HDM, sometimes called Type-Ill 2HDM 42|. In this case, the Feynman rules for 
the charged Higgs-boson couplings to quarks, with all particles incoming, are 



H+udj : -if-V^.,(A>„^Pi-A>,^,P^), 



V2Mw 
where the family-dependent couplings A\ ^ read 



(66) 



Al 



A, 



mr 



1 + Cn.d^T^JS 



(67) 



Under the assumption that no new sources of CP violation apart from the complex phase 
in the CKM matrix are present, the coefficients Au^d and eu4 are real. The case eu^ = 
corresponds to the NFC situation, in which the Yukawa matrices of both Higgs doublets are 
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aligned in flavor space. LEP measurements of Rh constrain \Au\ to values below 0.3 and 
0.5 (0.78 and 1.35) for mn = 100 and 400 GeV at la (2(t), when Aa = 0. For opposite 
(same) signs of Au and A^, the average of BABAR, Belle and CLEO measurements of 
BK{B — 7- Xg'-f) allow for one (two) region(s) of A^ for given values of A^ and rriH- For 
Au = 0.3 and rriH = 100 GeV, both A^ G [0; 1] and [16; 18] are allowed, while for Au = 0.3 



and rriH = 400 GeV, both Ad G [0; 2.5] and [50;56] are allowed at 2a J42|. Since general color- 
singlet Higgs-boson couplings and (theoretically possible) color-octet Higgs bosons induce 
different QCD corrections, we will not study these scenarios numerically. In the literature, 
one may also find Type-III 2HDMs where no flavor symmetry is imposed and FCNCs are 



avoided by other methods, e.g. by the small mass of first and second generation quarks |43 |. 
These models then allow for the couplings of charged Higgs bosons to bottom and charm 
quarks, which induces a phenomenology that is different from the one studied in this paper. 

C. Predictions for various 2HDMs 

The calculation presented in the previous sections was performed in a generic way, which 
makes it possible to use the result for various models with charged Higgs bosons. Out of the 
models mentioned in the last section, our calculation is in particular valid for the Type-I 
and Type-II 2HDMs. In this subsection, we perform a numerical analysis for a set of typical 
collider scenarios for these 2HDMs. As in these models the scattering matrix element is 
directly proportional to the Higgs-top-bottom quark coupling even at NLO, the type of the 
model has no influence on kinematic distributions apart from their normalization to the 
total cross section. 

We therefore concentrate here on the total cross sections and on the uncertainties both 
from the variation of the renormalization and factorization scales and from the parton dis- 
tribution functions. For a better comparison, we analyze the total cross sections and their 
uncertainties by choosing the same values for the mass of the charged Higgs boson and 
for tan/3 in all scenarios, i.e. rriH = 300 GeV and tan/3 = 10. All relevant values are 
summarized in Tab. llVi 

In all scenarios, both at the Tevatron and at the LHC, the next-to-leading order correction 
is substantial, ranging from 57% at the Tevatron in the Type-I 2HDM to 38% at the LHC in 
the same model. Apart from enhancing the total cross section, including the NLO correction 
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TABLE IV: Total cross sections (in pb) for different 2HDMs at the Tevatron and at the LHC at 
leading order (LO) and at next-to-leading order (NLO) including the scale and PDF uncertainties. 
All scenarios assume the same parameters for better comparison, i.e. mn = 300 GeV and tan/3 = 
10. 



Scenario 


LO 


Scale unc. 


NLO 


Scale unc. 


PDF error 


Tevatron 2HDM-I 
Tevatron 2HDM-II 


3.229.10-6 
1.303.10-5 


+1.306.10-6(40%) 
-0.901.10-6(28%) 
+0.524.10-5(40%) 
-0.365.10-5(28%) 


6.218.10-6 
2.506.10-5 


+ 1.388.10-6(22%) 
-1.201.10-6(19%) 
+0.565.10-5(23%) 
-0.484.10-5(19%) 


+4.448.10-5(72%) 
-2.362.10-5(38%) 
+ 1.792.10-5(72%) 
-0.952.10-5(38%) 


LHC 2HDM-I 
LHC 2HDM-n 


1.577.10-3 
6.366.10-3 


+0.379.10-3(24%) 
-0.304.10-3(19%) 
+1.514.10-3(24%) 
-1.237.10-3(19%) 


2.189.10-3 
8.821.10-3 


+0.162.10-3(7%) 
-0.199.10-3(9%) 
+0.651.10-3(7%) 
-0.802.10-3(9%) 


+0.356.10-3(16%) 
-0.304.10-3(14%) 
+ 1.433.10-3(16%) 
-1.223.10-3(14%) 



reduces the theoretical error defined as the scale uncertainty of the cross section. The 
scale uncertainty is obtained by varying both the renormalization and factorization scales 
simultaneously in the interval 

rrit + tuh 



< fi < rrit + rriH ■ 



(68) 



At leading order, the strong scale dependence conies from the strong coupling constant and 
from the Yukawa coupling in the tree-level amplitude. Including higher-order corrections, 
this uncertainty is dramatically reduced in some scenarios. 

Another large source of error stems from the parton distribution functions. We use 
the CTIO NLO PDF set with its error PDF sets to determine the error coming from the 
uncertainty contained in determining the parton content of the colliding hadrons. The 
process considered here is extremely sensitive to the gluon distribution function through 
having a gluon in the initial state and through having a heavy-quark initial state, which is 
radiatively generated from the gluon PDF. Moreover, the production of a heavy Higgs boson 
in association with a top quark probes the higher x content of the initial-state (anti-)proton. 
The values of Bjorken-x probed can be expressed as 



_ (fci + /C2)^ {mt + mn) 



(69) 

s s 

which at the Tevatron leads to typical values of a; ~ 0.3. This is exactly the region where 
the gluon PDF is poorly known, which translates into large PDF uncertainties on the cross 
section at the Tevatron. At the LHC, the Bjorken-x probed is x ~ 0.1, and the PDF 
uncertainties are therefore much smaller. 
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D. Checks of the NLO calculation and comparisons with POWHEG 

As a check of the numerical implementation of our analytical results, we have compared 
our complete NLO calculation obtained with the Catani-Seymour dipole formalism with 
the one performed previously with a phase-space slicing method using a single invariant- 
mass cutoff [12|, which had in turn been found to a gree with a calculation using a two 



(soft and collinear) cutoff phase-space slicing method ll|]. We found good agreement for 
all differential and total cross sections studied, but refrain from showing the corresponding 
figures here, since the fixed-order results are well-known. 

For the remainder of the analysis, we will constrain ourselves to the Type-II 2HDM, as 
the kinematic distributions have the same features in both Type-I and Type-II 2HDMs. In 
all of our discussion, we consider three collider scenarios: 

• Tevatron, y/S = 1.96 TeV, 

• LHC, v^ = 7 TeV, and 

• LHC, v^ = 14 TeV. 

Moreover, in the comparison of our NLO calculation with our implementation of its relevant 
parts in the POWHEG BOX, we assume rriH = 300 GeV and tan/3 = 10. The results are 
shown in Fig. IHfor the Tevatron with -\A5 = 1.96 TeV and Figs. [5] and |6] for the LHC with 
a center-of-mass energy of y/S = 7 and 14 TeV, respectively. 

If we concentrate first on the transverse-momentum {pr, left) and rapidity {y, right) dis- 
tributions of the charged Higgs boson (top) and top quark (center) individually, we observe 
good agreement in absolute normalization and shape for all three collider scenarios, indepen- 
dently if a parton shower is matched to the NLO calculation or not. This corresponds to the 
well-known fact that these distributions are largely insensitive to soft or collinear radiation, 
in particular from the initial state, and this can therefore be seen as a further consistency 
test of our calculations. Soft radiation becomes relevant in all three collider scenarios when 
we consider the azimuthal opening angle of the top-Higgs pair (bottom right), where the 
singularity occurring at NLO in back-to-back kinematics at A0 = vr is regularized and re- 
summed by the parton showers. This holds also for the p^-distribution of the top-Higgs pair 
(bottom left), which diverges perturbatively at p^ = GeV and even turns negative at the 
LHC. 
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POWHEG - Tevatron, Vs = L96TeV, 2HDM-IL mH = 300 GeV, tan6 = 10 
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FIG. 4: Distributions in transverse momentum -pT (top left) and rapidity y (top right) of the 

charged Higgs boson, pT (center left) and y (center right) of the top quark, as well as "pT (bottom 

left) and azimuthal opening angle A</> (bottom right) of the tH^ system produced at the Tevatron 

with \fS = 1.96 TeV. We compare the NLO predictions without (blue) and with matching to the 

PYTHIA (black) and HERWIG (red) parton showers using POWHEG in the Type-II 2HDM with 

tan/3 = 10 and niH = 300 GeV. 
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POWHEG - LHC, Vs = 7TeV, 2HDM-II, mH = 300 GeV, tan6 = 10 
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FIG. 5: Same as Fig. Hat the LHC with VS = 7 TeV. 

An advantage of the POWHEG method is that it can also provide events including first 
radiation only in the form of an event file according to the Les Houches format (LHEF), 
making them independent of the parton shower. We therefore compare in Figs. U\M the 
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POWHEG - LHC, Vs = 14TeV, 2HDM-IL mH = 300 GeV, tan6 = 10 
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FIG. 6: Same as Fig. Hat the LHC with VS = U TeV. 

distributions obtained from these files to those obtained with NLO accuracy for the same 
set of parameters as in Figs. HHEl As one can clearly see, they lie within the NLO scale 
uncertainty band, showing that the difference comes from terms beyond NLO accuracy. 
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POWHEG - Tevatron, Vs = L96TeV, 2HDM-IL mH = 300 GeV, tan6 = 10 
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FIG. 7: Distributions in transverse momentum pT (top left) and rapidity y (top right) of the 
charged Higgs boson, pT (center left) and y (center right) of the top quark, as well as pT (bottom 
left) and azimuthal opening angle A</> (bottom right) of the tH^ system produced at the Tevatron 
with v5^ = 1.96 TeV. We compare the NLO scale uncertainty band (blue) the POWHEG result 
including first radiation only (red). 
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POWHEG - LHC, Vs = 7TeV, 2HDM-II, mH = 300 GeV, tan6 = 10 



6.0-10" 



> 

CD 
CD 



L4.0-10" 



to2.0-10" 




BLUEs NLO 
RED: POWHEG 



_L 




100 200 300 

Pt H (GeV) 



6.0-10" 



> 

ID 



L4.0-10" 



ID2.0-10" 



100 200 300 

Pj top (GeV) 



> 

CD 



100 200 300 

Pt tH (GeV) 



b 



400 




400 




400 




A0 tH 



FIG. 8: Same as Fig. [7] at the LHC with VS = 7 TeV. 



This provides a good consistency check of the matching procedure. 
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POWHEG - LHC, Vs = 14TeV, 2HDM-IL mH = 300 GeV, tan6 = 10 
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FIG. 9: Same as Fig. Hat the LHC with ^5 = 14 TeV. 
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E. POWHEG predictions with HERWIG and PYTHIA parton showers 

In Figs. HHHl we also show two different predictions with POWHEG coupled either to 
the angularly-ordered HERWIG or to the virtuality-ordered PYTHIA parton shower. The 
agreement of the HERWIG and PYTHIA results is in general very good. They differ only 
slightly in the px distributions of the top-Higgs pair, where the PYTHIA pT-distribution is 
a little bit harder, in particular at the Tevatron. 

F. POWHEG comparison with MATCHIG 

As described in Sec. Ull the production of charged Higgs bosons and top quarks proceeds 
at LO through the process bg — > H^t, while at NLO the process gg — > H^tb appears. 
The latter implies the creation of a virtual initial 6-quark, which may either occur in the 
perturbative part of the calculation or is resummed into a 6-quark PDF. In the full NLO 
calculation, the separation is achieved through the factorization procedure and induces a 
dependence on the factorization scale fip- 

Before schemes to match parton showers with full NLO calculations were developed, the 
importance of the contribution of this particular two-to-three_process and the perturbative 



origin of the 6-quark density had already been recognized (25|. It had been proposed to 
supplement the LO calculation by this particular two-to-three process and to remove the 
overlap by subtracting the doubly counted (DC) term 

O-DC = / dXafb'iXa^lJ-F) dXb fg{Xb, IJ,F)a^'^ {pi, P2) + {Xa <-^ Xb), (70) 

Jo Jo 

where fb'{x,fip) is the LO 6-quark density given by 

Mx,^l) - f^ln4 l-P.9i^)h (-,/^l) (71) 

with Pqg{z) the g ^f q splitting function, fg{x,fip) the gluon PDF, and z the longitudinal 
gluon momentum fraction taken by the 6-quark. The two-to-three and double- counting 
processes had been implemented in an addition to PYTHIA called MATCHIG. 

With our full NLO calculation matched to PYTHIA within the POWHEG BOX, it is now 
possible to compare the two approaches numerically. The results are shown in Fig. [101 Since 
the normalization of the MATCHIG prediction is still effectively of LO, we have normalized 
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MATCHIG-POWHEG comparison - LHC, Vs = 7TeV, 2HDM-II, mH = 300 GeV, tanB^lO 
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FIG. 10: Distributions in transverse momentum pT (top left) and rapidity y (top right) of the 

charged Higgs boson, pT (center left) and y (center right) of the top quark, as well as pT (bottom 

left) and azimuthal opening angle A(j) (bottom right) of the tH^ system produced at the LHC with 

yo = 7 TeV. We compare the tree-level predictions matched to PYTHIA using MATCHIG (black) 

with our NLO calculation matched to PYTHIA (red) and HERWIG (blue) using POWHEG. Ah 

distributions have been normalized to the respective total cross sections. 
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all distributions to their respective total cross sections in order to emphasize the shapes of 
the distributions. One observes that when both the MATCHIG (black) and POWHEG (red) 
predictions are matched to the PYTHIA parton shower, there is very little difference, even 
at low Pt and large A0 of the top-Higgs pair. Only at large px and small A0 the differences 
become sizable, which can be attributed to the fact that MATCHIG includes only one of 
the four classes of real-emission processes, while our POWHEG prediction includes also the 
quark-initiated real-emission processes. Let us emphasize again that while the spectra are 
already quite well described with MATCHIG, their normalization is only accurate to LO 
and not NLO as in POWHEG. 

G. Comparison with MC@NLO 

In a recent publication, two of us and a number of other authors have matched a NLO 
calculation performed with the FKS subtraction formalism to the HERWIG PS with the 



MC@NLO method [20|. It is therefore mandatory that we compare in this paper this 
previous work with our new POWHEG implementation, which we do in Fig. [11] Note that 
here we employ a value of tan /3 = 30 as in the MC@NLO publication. In both calculations, 
we use the HERWIG PS in order to emphasize possible differences in the matching methods 
and not those in the parton shower. We also normalize the differential cross sections again 
to the total cross section for a better comparison of the shapes of the distributions. 

As in the other comparisons, the rapidity distributions of the charged Higgs boson (top 
right) and the top quark (center right) show little variation, confirming the consistency 
of the two calculations. However, the corresponding p^-spectra (top and center left) are 
slightly harder with the MC@NLO matching than in POWHEG. This behaviour is known 
from other processes 
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44| . It is less pronounced in the pr-distribution of the top-Higgs 
pair, shown on a logarithmic scale (bottom left). Since we are using the HERWIG PS, the 
rise at small azimuthal angle A0 (bottom right) is not very strong with MC@NLO and only 
slightly more so with POWHEG. In total, all of these differences are similarly small in the 
production of a top quark with a IV-boson ^] and with a charged Higgs boson at the LHC. 
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MC@NLO-POWHEG comparison - LHC, Vs = 14TeV, 2HDM-II, mH = 300 GeV, tan6 = 30 
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FIG. 11: Distributions in transverse momentum pT (top left) and rapidity y (top right) of the 

charged Higgs boson, px (center left) and y (center right) of the top quark, as well as pT (bottom 

left) and azimuthal opening angle A(j) (bottom right) of the tH^ system produced at the LHC 

with VJ = 14 TeV. We compare the NLO predictions with matching to the HERWIG parton 

showers using POWHEG (red) and MC@NLO (black) in the Type-II 2HDM with tan/3 = 30 and 

ruH = 300 GeV. All distributions have been normalized to the respective total cross sections. 
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H. Diagram Removal, Diagram Subtraction, and no subtraction 

If the charged Higgs boson was hghter than the top quark, it would dominantly be created 
in top-pair production and the decay of an (anti-)top quark into it. As discussed above, 
one must then find a suitable definition to separate this process from the associated top- 
Higgs production discussed in this paper. In addition to the Diagram Removal (DR) and 
Diagram Subtraction (DS) methods discussed above, we introduce here also the option of 
not removing or subtracting anything from the associated production, but simply retaining 
the total production cross section, which then allows for the removal of fully simulated events 
near the resonance region and replacing them with events obtained, e.g., with a full NLO 
implementation of ti production. The results are shown in Fig. [T21 

The rapidity distributions of the charged Higgs boson (top right) and top quark (center 
right) show again little sensitivity to the different theoretical approaches. However, the 
PT-distribution of the charged Higgs boson (top left) is somewhat softer and the one of the 
top quark (center left) considerably harder without removal or subtraction, as the difference 
describes the distributions of the lighter decay product and the heavier decaying particle, 
respectively. The p^-distribution of the top-Higgs pair (bottom left) is significantly harder 
(note again the logarithmic scale) and its maximum moves from px = 20 to 70 GeV, indi- 
cating that the transverse momentum of the pair is balanced by a hard object, i.e. the fast 
additional 6-quark jet, in the other hemisphere. This also allows the top-Higgs pair to move 
closer together in azimuthal angle (bottom right). 

The theoretical pros and cons and the numerical differences of Diagram Removal and 



Diagram Subtraction have been discussed extensively above and also elsewhere 20|. It is 
clear from Fig. [12] that the numerical difference of DR vs. DS is much less pronounced than 
the difference of both with respect to no removal or subtraction at all. We emphasize that 
the total cross section is continuous across the rriH = Trit threshold in all three schemes (see 



also Ref. |27|) 



The differences of POWHEG and MC@NLO are small for itih < frit in both the DR and 
DS schemes, as can be seen when comparing Figs. [12] and [TH This coincides nicely with our 
observation above that these differences should be as small as in the associated production 
of VT-bosons and top quarks 
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DIAGRAM REMOVAL vs SUBTRACTION - LHC, Vs = 14TeV, 2HDM-II, mH = 100 GeV, tan6 = 30 



0.010 



BLUE: ail 
RED: POWHEG-DR 
BLACK: POWHEG-DS 




200 300 

Pt H (6eV) 



100 200 300 

Pj top (GgV) 



0.05000 

> 0.01000 
^ 0.00500 

:§" 0.00100 

-g 0.00050 

b 

0.00010 
0.00005 



400 





0.008 - 


— 


> 

CD 


0.006 - 


A - 


" 


- 


r ""n 


Q. 

■o 


0.004 - 


/ \\ 


b 


0.002 - 


r i>Su ~ 




: 


Ti"h_ 




= 


^~'~'=b^^!l*~M_ 




n nnn 


. . , . 1 . . , . 1 . , 7r"^^=^=^= 



400 




100 200 300 

Pt tH (GeV) 



400 




FIG. 12: Distributions in transverse momentum pT (top left) and rapidity y (top right) of the 
charged Higgs boson, pT (center left) and y (center right) of the top quark, as well as pT (bottom 
left) and azimuthal opening angle A(j) (bottom right) of the tH~ system produced at the LHC 
with \/S = 14 TeV. We compare the NLO predictions matched to the HERWIG parton shower 
using POWHEG with Diagram Removal (red). Diagram Subtraction (black), and without removing 
or subtracting anything (blue) in the Type-II 2HDM with tan/3 = 30 and itlh = 100 GeV. All 
distributions have been normalized to the respective total cross sections. 
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DIAGRAM REMOVAL POWHEG vs MC@NLO - LHC, Vs = 14TeV, 2HDM-II, mH = 100 GeV, tan6 = 30 
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FIG. 13: Same as Fig. [TTl but for a light charged Higgs boson of mass mn = 100 GeV and using 
the DR method. 

V. CONCLUSION 

In this paper, we presented a new NLO calculation of the associated production of charged 
Higgs bosons and top quarks at hadron colliders using the Catani-Seymour dipole subtraction 
formalism and matched it to parton showers with the POWHEG method. We discussed the 
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DIAGRAM SUBTRACTION POWHEG vs MC@NLO - LHC, Vs = 14TeV, 2HDM-II, mH=100 GeV, tan6=30 
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FIG. 14: Same as Fig. [TTl but for a light charged Higgs boson of mass mn = 100 GeV and using 
the DS method. 

different types of 2HDMs as well as the corresponding current experimental constraints and 
provided, for specific benchmark values of the charged Higgs-boson mass and the ratio of 
the two Higgs VEVs tan/3, the central values, scale, and PDF uncertainties of the total 
cross sections at the Tevatron and LHC in tabular form for future reference. As expected, 
the scale uncertainty was considerably reduced from up to ±100% at LO to less than ±15% 
at NLO. However, the PDF uncertainty, estimated with the CTIO set of global analyses. 
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remained quite substantial, in particular at the Tevatron, where high momentum fractions 
of the gluons and 6-quarks in the protons and antiprotons are probed. 

For the differential cross sections, we established good numerical agreement of our full 
NLO calculation with previous calculations. We then performed detailed comparisons of 
our new POWHEG implementation with the purely perturbative result, with PYTHIA or 
HERWIG parton showers, with a LO calculation matched to the PYTHIA parton shower 
using MATCHIG, and with a NLO calculation matched to the HERWIG parton shower 
using MC@NLO. 

While the transverse-momentum distributions and the relatively central rapidity distri- 
butions of the charged Higgs boson and top quark individually showed little sensitivity 
to the existence and type of parton showers, the transverse-momentum distribution of the 
top-Higgs pair depended quite strongly on the different theoretical approaches as expected. 
This was also true for the distribution in the azimuthal angle of the top-Higgs pair. For 
scenarios in which the charged Higgs boson is lighter than the top quark, we implemented in 
POWHEG in addition to the previously proposed Diagram Removal and Diagram Subtrac- 
tion schemes the possibility to retain the full cross section and replace the simulated events 
in the resonance region with a full NLO Monte Carlo for top-quark pair production. 

It will now be very interesting to observe the impact of our work on the experimental 
search for charged Higgs bosons. The numerical code and technical support is, of course, 
available from the authors. 
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